home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / PRISM.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-18  |  42.2 KB  |  1,949 lines
  1. /****************************************************************************
  2. *                   prism.c
  3. *
  4. *  This module implements functions that manipulate prisms.
  5. *
  6. *  This module was written by Dieter Bayer [DB].
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file. If
  16. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  18. *  Forum.  The latest version of POV-Ray may be found there as well.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. /****************************************************************************
  27. *
  28. *  Explanation:
  29. *
  30. *    The prism primitive is defined by a set of points in 2d space which
  31. *    are interpolated by linear, quadratic, or cubic splines. The resulting
  32. *    2d curve is swept along a straight line to form the final prism object.
  33. *
  34. *    All calculations are done in the object's (u,v,w)-coordinate system
  35. *    with the (w)-axis being the sweep axis.
  36. *
  37. *    One spline segment in the (u,v)-plane is given by the equations
  38. *
  39. *      fu(t) = Au * t * t * t + Bu * t * t + Cu * t + Du  and
  40. *      fv(t) = Av * t * t * t + Bv * t * t + Cv * t + Dv,
  41. *
  42. *    with the parameter t ranging from 0 to 1.
  43. *
  44. *    To intersect a ray R = P + k * D transformed into the object's
  45. *    coordinate system with the linear swept prism object, the equation
  46. *
  47. *      Dv * fu(t) - Du * fv(t) - (Pu * Dv - Pv * Du) = 0
  48. *
  49. *    has to be solved for t. For valid intersections (0 <= t <= 1)
  50. *    the corresponding k can be calculated from equation
  51. *
  52. *      Pu + k * Du = fu(t) or Pv + k * Dv = fv(t).
  53. *
  54. *    In the case of conic sweeping the equation
  55. *
  56. *      (Pv * Dw - Pw * Dv) * fu(t) - (Pu * Dw - Pw * Du) * fv(t)
  57. *                                  + (Pu * Dv - Pv * Du) = 0
  58. *
  59. *    has to be solved for t. For valid intersections (0 <= t <= 1)
  60. *    the corresponding k can be calculated from equation
  61. *
  62. *      Pu + k * Du = (Pw + k * Dw) * fu(t) or
  63. *      Pv + k * Dv = (Pw + k * Dw) * fv(t).
  64. *
  65. *    Note that the polynomal to solve has the same degree as the
  66. *    spline segments used.
  67. *
  68. *    Note that Pu, Pv, Pw and Du, Dv, Dw denote the coordinates
  69. *    of the vectors P and D.
  70. *
  71. *  Syntax:
  72. *
  73. *    prism {
  74. *      [ linear_sweep | cubic_sweep ]
  75. *      [ linear_spline | quadratic_spline | cubic_spline ]
  76. *
  77. *      height1, height2,
  78. *      number_of_points
  79. *
  80. *      <P(0)>, <P(1)>, ..., <P(n-1)>
  81. *
  82. *      [ open ]
  83. *      [ sturm ]
  84. *    }
  85. *
  86. *    Note that the P(i) are 2d vectors.
  87. *
  88. *  ---
  89. *
  90. *  Ideas for the prism was taken from:
  91. *
  92. *    James T. Kajiya, "New techniques for ray tracing procedurally
  93. *    defined objects", Computer Graphics, 17(3), July 1983, pp. 91-102
  94. *
  95. *    Kirk ...
  96. *
  97. *  ---
  98. *
  99. *  May 1994 : Creation. [DB]
  100. *
  101. *****************************************************************************/
  102.  
  103. #include "frame.h"
  104. #include "povray.h"
  105. #include "vector.h"
  106. #include "povproto.h"
  107. #include "bbox.h"
  108. #include "matrices.h"
  109. #include "objects.h"
  110. #include "polysolv.h"
  111. #include "prism.h"
  112.  
  113.  
  114.  
  115. /*****************************************************************************
  116. * Local preprocessor defines
  117. ******************************************************************************/
  118.  
  119. /* Minimal intersection depth for a valid intersection. */
  120.  
  121. #define DEPTH_TOLERANCE 1.0e-4
  122.  
  123. /* |Coefficients| < COEFF_LIMIT are regarded to be 0. */
  124.  
  125. /*#define COEFF_LIMIT 1.0e-20 */
  126.  
  127. #define COEFF_LIMIT 1.0e-16 /*changed by CEY 11/18/95 */
  128.  
  129. /* Part of the prism hit. */
  130.  
  131. #define BASE_HIT   1
  132. #define CAP_HIT    2
  133. #define SPLINE_HIT 3
  134.  
  135.  
  136.  
  137. /*****************************************************************************
  138. * Local typedefs
  139. ******************************************************************************/
  140.  
  141. typedef struct Prism_Intersection_Structure PRISM_INT;
  142.  
  143. struct Prism_Intersection_Structure
  144. {
  145.   DBL d;  /* Distance of intersection point                  */
  146.   DBL w;  /* Paramter of intersection point on n-th spline   */
  147.   int n;  /* Number of segment hit                           */
  148.   int t;  /* Type of intersection: base/cap plane or segment */
  149. };
  150.  
  151.  
  152.  
  153. /*****************************************************************************
  154. * Static functions
  155. ******************************************************************************/
  156.  
  157. static int intersect_prism PARAMS((RAY *Ray, PRISM *Prism, PRISM_INT *Intersection));
  158. static int in_curve PARAMS((PRISM *Prism, DBL u, DBL v));
  159. static int test_rectangle PARAMS((VECTOR P, VECTOR D, DBL x1, DBL y1, DBL x2, DBL y2));
  160. static int   All_Prism_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  161. static int   Inside_Prism PARAMS((VECTOR point, OBJECT *Object));
  162. static void  Prism_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  163. static void  *Copy_Prism PARAMS((OBJECT *Object));
  164. static void  Translate_Prism PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  165. static void  Rotate_Prism PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  166. static void  Scale_Prism PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  167. static void  Transform_Prism PARAMS((OBJECT *Object, TRANSFORM *Trans));
  168. static void  Invert_Prism PARAMS((OBJECT *Object));
  169. static void  Destroy_Prism PARAMS((OBJECT *Object));
  170.  
  171.  
  172. /*****************************************************************************
  173. * Local variables
  174. ******************************************************************************/
  175.  
  176. static METHODS Prism_Methods =
  177. {
  178.   All_Prism_Intersections,
  179.   Inside_Prism, Prism_Normal,
  180.   Copy_Prism,
  181.   Translate_Prism, Rotate_Prism,
  182.   Scale_Prism, Transform_Prism, Invert_Prism, Destroy_Prism
  183. };
  184.  
  185.  
  186.  
  187. /*****************************************************************************
  188. *
  189. * FUNCTION
  190. *
  191. *   All_Prism_Intersections
  192. *
  193. * INPUT
  194. *
  195. *   Object      - Object
  196. *   Ray         - Ray
  197. *   Depth_Stack - Intersection stack
  198. *   
  199. * OUTPUT
  200. *
  201. *   Depth_Stack
  202. *   
  203. * RETURNS
  204. *
  205. *   int - TRUE, if a intersection was found
  206. *   
  207. * AUTHOR
  208. *
  209. *   Dieter Bayer
  210. *   
  211. * DESCRIPTION
  212. *
  213. *   Determine ray/prism intersection and clip intersection found.
  214. *
  215. * CHANGES
  216. *
  217. *   May 1994 : Creation.
  218. *
  219. ******************************************************************************/
  220.  
  221. static int All_Prism_Intersections(Object, Ray, Depth_Stack)
  222. OBJECT *Object;
  223. RAY *Ray;
  224. ISTACK *Depth_Stack;
  225. {
  226.   int i, max_i, Found;
  227.   PRISM_INT *I;
  228.   VECTOR IPoint;
  229.  
  230.   Found = FALSE;
  231.  
  232.   /* Allocate memory for intersections. */
  233.  
  234.   I = (PRISM_INT *)POV_MALLOC((((PRISM *)Object)->Number+2)*sizeof(PRISM_INT), "prism intersection array");
  235.  
  236.   max_i = intersect_prism(Ray, (PRISM *)Object, &I[0]);
  237.  
  238.   if (max_i)
  239.   {
  240.     for (i = 0; i < max_i; i++)
  241.     {
  242.       if ((I[i].d > DEPTH_TOLERANCE) && (I[i].d < Max_Distance))
  243.       {
  244.         VEvaluateRay(IPoint, Ray->Initial, I[i].d, Ray->Direction);
  245.  
  246.         if (Point_In_Clip(IPoint, Object->Clip))
  247.         {
  248.           push_entry_i1_i2_d1(I[i].d,IPoint,Object,I[i].t,I[i].n,I[i].w,Depth_Stack);
  249.  
  250.           Found = TRUE;
  251.         }
  252.       }
  253.     }
  254.   }
  255.  
  256.   POV_FREE (I);
  257.  
  258.   return(Found);
  259. }
  260.  
  261.  
  262.  
  263. /*****************************************************************************
  264. *
  265. * FUNCTION
  266. *
  267. *   intersect_prism
  268. *
  269. * INPUT
  270. *
  271. *   Ray   - Ray
  272. *   Prism - Prism
  273. *   Int   - Prism intersection structure
  274. *   
  275. * OUTPUT
  276. *
  277. *   Int
  278. *   
  279. * RETURNS
  280. *
  281. *   int - Number of intersections found
  282. *   
  283. * AUTHOR
  284. *
  285. *   Dieter Bayer
  286. *   
  287. * DESCRIPTION
  288. *
  289. *   Determine ray/prism intersection.
  290. *
  291. *   NOTE: Order reduction cannot be used.
  292. *
  293. * CHANGES
  294. *
  295. *   May 1994 : Creation.
  296. *
  297. ******************************************************************************/
  298.  
  299. static int intersect_prism(Ray, Prism, Intersection)
  300. RAY *Ray;
  301. PRISM *Prism;
  302. PRISM_INT *Intersection;
  303. {
  304.   int i, j, n;
  305.   DBL k, u, v, w, h, len;
  306.   DBL x[4], y[3];
  307.   DBL k1, k2, k3;
  308.   VECTOR P, D;
  309.   PRISM_SPLINE_ENTRY Entry;
  310.   DBL rdiv;
  311.   
  312.   /* Don't test degenerate prisms. */
  313.  
  314.   if (Test_Flag(Prism, DEGENERATE_FLAG))
  315.   {
  316.     return(FALSE);
  317.   }
  318.  
  319.   Increase_Counter(stats[Ray_Prism_Tests]);
  320.  
  321.   /* Init intersection structure. */
  322.  
  323.   Intersection->d =
  324.   Intersection->w = 0.0;
  325.   Intersection->n =
  326.   Intersection->t = 0;
  327.  
  328.   /* Transform the ray into the prism space */
  329.  
  330.   MInvTransPoint(P, Ray->Initial, Prism->Trans);
  331.  
  332.   MInvTransDirection(D, Ray->Direction, Prism->Trans);
  333.  
  334.   VLength(len, D);
  335.  
  336.   VInverseScaleEq(D, len);
  337.  
  338.   /* Test overall bounding rectangle. */
  339.  
  340. #ifdef PRISM_EXTRA_STATS
  341.   Increase_Counter(stats[Prism_Bound_Tests]);
  342. #endif
  343.  
  344.   if (((D[X] >= 0.0) && (P[X] > Prism->x2)) ||
  345.       ((D[X] <= 0.0) && (P[X] < Prism->x1)) ||
  346.       ((D[Z] >= 0.0) && (P[Z] > Prism->y2)) ||
  347.       ((D[Z] <= 0.0) && (P[Z] < Prism->y1)))
  348.   {
  349.     return(FALSE);
  350.   }
  351.  
  352. #ifdef PRISM_EXTRA_STATS
  353.   Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  354. #endif
  355.  
  356.   /* Number of intersections found. */
  357.  
  358.   i = 0;
  359.  
  360.   /* What kind of sweep is used? */
  361.  
  362.   switch (Prism->Sweep_Type)
  363.   {
  364.     /*************************************************************************
  365.     * Linear sweep.
  366.     **************************************************************************/
  367.  
  368.     case LINEAR_SWEEP :
  369.  
  370.       rdiv = (1.0 / D[Y]);
  371.  
  372.       if (fabs(D[Y]) < EPSILON)
  373.       {
  374.         if ((P[Y] < Prism->Height1) || (P[Y] > Prism->Height2))
  375.         {
  376.           return(FALSE);
  377.         }
  378.       }
  379.       else
  380.       {
  381.         if (Test_Flag(Prism, CLOSED_FLAG))
  382.         {
  383.           /* Intersect ray with the cap-plane. */
  384. /*
  385.           k = (Prism->Height2 - P[Y]) / D[Y];
  386. */ 
  387.           k = (Prism->Height2 - P[Y]) * rdiv;
  388.  
  389.           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  390.           {
  391.             u = P[X] + k * D[X];
  392.             v = P[Z] + k * D[Z];
  393.  
  394.             if (in_curve(Prism, u, v))
  395.             {
  396.               Intersection[i].t   = CAP_HIT;
  397.               Intersection[i++].d = k / len;
  398.             }
  399.           }
  400.  
  401.           /* Intersect ray with the base-plane. */
  402. /*
  403.           k = (Prism->Height1 - P[Y]) / D[Y];
  404. */
  405.           k = (Prism->Height1 - P[Y]) * rdiv;
  406.  
  407.           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  408.           {
  409.             u = P[X] + k * D[X];
  410.             v = P[Z] + k * D[Z];
  411.  
  412.             if (in_curve(Prism, u, v))
  413.             {
  414.               Intersection[i].t   = BASE_HIT;
  415.               Intersection[i++].d = k / len;
  416.             }
  417.           }
  418.         }
  419.       }
  420.  
  421.       /* Intersect ray with all spline segments. */
  422.  
  423.       if ((fabs(D[X]) > EPSILON) || (fabs(D[Z]) > EPSILON))
  424.       {
  425.         for (j = 0; j < Prism->Number; j++)
  426.         {
  427.           Entry = Prism->Spline->Entry[j];
  428.  
  429.           /* Test spline's bounding rectangle (modified Cohen-Sutherland). */
  430.  
  431. #ifdef PRISM_EXTRA_STATS
  432.           Increase_Counter(stats[Prism_Bound_Tests]);
  433. #endif
  434.  
  435.           if (((D[X] >= 0.0) && (P[X] > Entry.x2)) ||
  436.               ((D[X] <= 0.0) && (P[X] < Entry.x1)) ||
  437.               ((D[Z] >= 0.0) && (P[Z] > Entry.y2)) ||
  438.               ((D[Z] <= 0.0) && (P[Z] < Entry.y1)))
  439.           {
  440.             continue;
  441.           }
  442.  
  443.           /* Number of roots found. */
  444.  
  445.           n = 0;
  446.  
  447.           switch (Prism->Spline_Type)
  448.           {
  449.             case LINEAR_SPLINE :
  450.  
  451. #ifdef PRISM_EXTRA_STATS
  452.               Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  453. #endif
  454.  
  455.               /* Solve linear equation. */
  456.  
  457.               x[0] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
  458.  
  459.               x[1] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
  460.  
  461.               if (fabs(x[0]) > EPSILON)
  462.               {
  463.                 y[n++] = -x[1] / x[0];
  464.               }
  465.  
  466.               break;
  467.  
  468.             case QUADRATIC_SPLINE :
  469.  
  470. #ifdef PRISM_EXTRA_STATS
  471.               Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  472. #endif
  473.  
  474.               /* Solve quadratic equation. */
  475.  
  476.               x[0] = Entry.B[X] * D[Z] - Entry.B[Y] * D[X];
  477.  
  478.               x[1] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
  479.  
  480.               x[2] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
  481.  
  482.               n = Solve_Polynomial(2, x, y, FALSE, 0.0);
  483.  
  484.               break;
  485.  
  486.             case CUBIC_SPLINE :
  487.  
  488.               if (test_rectangle(P, D, Entry.x1, Entry.y1, Entry.x2, Entry.y2))
  489.               {
  490. #ifdef PRISM_EXTRA_STATS
  491.                 Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  492. #endif
  493.  
  494.                 /* Solve cubic equation. */
  495.  
  496.                 x[0] = Entry.A[X] * D[Z] - Entry.A[Y] * D[X];
  497.  
  498.                 x[1] = Entry.B[X] * D[Z] - Entry.B[Y] * D[X];
  499.  
  500.                 x[2] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
  501.  
  502.                 x[3] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
  503.  
  504.                 n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
  505.               }
  506.  
  507.               break;
  508.           }
  509.  
  510.           /* Test roots for valid intersections. */
  511.  
  512.           while (n--)
  513.           {
  514.             w = y[n];
  515.  
  516.             if ((w >= 0.0) && (w <= 1.0))
  517.             {
  518.               if (fabs(D[X]) > EPSILON)
  519.               {
  520.                 k = (w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X] - P[X]) / D[X];
  521.               }
  522.               else
  523.               {
  524.                 k = (w * (w * (w * Entry.A[Y] + Entry.B[Y]) + Entry.C[Y]) + Entry.D[Y] - P[Z]) / D[Z];
  525.               }
  526.  
  527.               /* Verify that intersection height is valid. */
  528.  
  529.               h = P[Y] + k * D[Y];
  530.  
  531.               if ((h >= Prism->Height1) && (h <= Prism->Height2))
  532.               {
  533.                 Intersection[i].t   = SPLINE_HIT;
  534.                 Intersection[i].n   = j;
  535.                 Intersection[i].w   = w;
  536.                 Intersection[i++].d = k / len;
  537.               }
  538.             }
  539.           }
  540.         }
  541.       }
  542.  
  543.       break;
  544.  
  545.     /*************************************************************************
  546.     * Conic sweep.
  547.     **************************************************************************/
  548.  
  549.     case CONIC_SWEEP :
  550.  
  551.       if (fabs(D[Y]) < EPSILON)
  552.       {
  553.         if ((P[Y] < Prism->Height1) || (P[Y] > Prism->Height2))
  554.         {
  555.           return(FALSE);
  556.         }
  557.       }
  558.       else
  559.       {
  560.         /* Intersect ray with the cap-plane. */
  561.  
  562.         k = (Prism->Height2 - P[Y]) / D[Y];
  563.  
  564.         if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  565.         {
  566.  
  567.           rdiv = (1.0 / Prism->Height2);
  568.           u = (P[X] + k * D[X]) * rdiv;
  569.           v = (P[Z] + k * D[Z]) * rdiv;
  570. /*
  571.           u = (P[X] + k * D[X]) / Prism->Height2;
  572.           v = (P[Z] + k * D[Z]) / Prism->Height2;
  573. */
  574.           if (in_curve(Prism, u, v))
  575.           {
  576.             Intersection[i].t   = CAP_HIT;
  577.             Intersection[i++].d = k / len;
  578.           }
  579.         }
  580.  
  581.         /* Intersect ray with the base-plane. */
  582.  
  583.         if (Prism->Height1 > 0.0)
  584.         {
  585.           k = (Prism->Height1 - P[Y]) / D[Y];
  586.  
  587.           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  588.           {
  589.  
  590.             rdiv = (1.0 / Prism->Height1);
  591.             u = (P[X] + k * D[X]) * rdiv;
  592.             v = (P[Z] + k * D[Z]) * rdiv;
  593. /*
  594.             u = (P[X] + k * D[X]) / Prism->Height1;
  595.             v = (P[Z] + k * D[Z]) / Prism->Height1;
  596. */
  597.             if (in_curve(Prism, u, v))
  598.             {
  599.               Intersection[i].t   = BASE_HIT;
  600.               Intersection[i++].d = k / len;
  601.             }
  602.           }
  603.         }
  604.       }
  605.  
  606.       /* Precompute ray-only dependant constants. */
  607.  
  608.       k1 = P[Z] * D[Y] - P[Y] * D[Z];
  609.  
  610.       k2 = P[Y] * D[X] - P[X] * D[Y];
  611.  
  612.       k3 = P[X] * D[Z] - P[Z] * D[X];
  613.  
  614.       /* Intersect ray with the spline segments. */
  615.  
  616.       if ((fabs(D[X]) > EPSILON) || (fabs(D[Z]) > EPSILON))
  617.       {
  618.         for (j = 0; j < Prism->Number; j++)
  619.         {
  620.           Entry = Prism->Spline->Entry[j];
  621.  
  622.           /* Test spline's bounding rectangle (modified Cohen-Sutherland). */
  623.  
  624.           if (((D[X] >= 0.0) && (P[X] > Entry.x2)) ||
  625.               ((D[X] <= 0.0) && (P[X] < Entry.x1)) ||
  626.               ((D[Z] >= 0.0) && (P[Z] > Entry.y2)) ||
  627.               ((D[Z] <= 0.0) && (P[Z] < Entry.y1)))
  628.           {
  629.             continue;
  630.           }
  631.  
  632.           /* Number of roots found. */
  633.  
  634.           n = 0;
  635.  
  636.           switch (Prism->Spline_Type)
  637.           {
  638.             case LINEAR_SPLINE :
  639.  
  640.               /* Solve linear equation. */
  641.  
  642.               x[0] = Entry.C[X] * k1 + Entry.C[Y] * k2;
  643.  
  644.               x[1] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
  645.  
  646.               if (fabs(x[0]) > EPSILON)
  647.               {
  648.                 y[n++] = -x[1] / x[0];
  649.               }
  650.  
  651.               break;
  652.  
  653.             case QUADRATIC_SPLINE :
  654.  
  655.               /* Solve quadratic equation. */
  656.  
  657.               x[0] = Entry.B[X] * k1 + Entry.B[Y] * k2;
  658.  
  659.               x[1] = Entry.C[X] * k1 + Entry.C[Y] * k2;
  660.  
  661.               x[2] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
  662.  
  663.               n = Solve_Polynomial(2, x, y, FALSE, 0.0);
  664.  
  665.               break;
  666.  
  667.             case CUBIC_SPLINE :
  668.  
  669.               /* Solve cubic equation. */
  670.  
  671.               x[0] = Entry.A[X] * k1 + Entry.A[Y] * k2;
  672.  
  673.               x[1] = Entry.B[X] * k1 + Entry.B[Y] * k2;
  674.  
  675.               x[2] = Entry.C[X] * k1 + Entry.C[Y] * k2;
  676.  
  677.               x[3] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
  678.  
  679.               n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
  680.  
  681.               break;
  682.           }
  683.  
  684.           /* Test roots for valid intersections. */
  685.  
  686.           while (n--)
  687.           {
  688.             w = y[n];
  689.  
  690.             if ((w >= 0.0) && (w <= 1.0))
  691.             {
  692.               k = w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X];
  693.  
  694.               h = D[X] - k * D[Y];
  695.  
  696.               if (fabs(h) > EPSILON)
  697.               {
  698.                 k = (k * P[Y] - P[X]) / h;
  699.               }
  700.               else
  701.               {
  702.                 k = w * (w * (w * Entry.A[Y] + Entry.B[Y]) + Entry.C[Y]) + Entry.D[Y];
  703.  
  704.                 h = D[Z] - k * D[Y];
  705.  
  706.                 if (fabs(h) > EPSILON)
  707.                 {
  708.                   k = (k * P[Y] - P[Z]) / h;
  709.                 }
  710.                 else
  711.                 {
  712.                   /* This should never happen! */
  713.                   continue;
  714.                 }
  715.               }
  716.  
  717.               /* Verify that intersection height is valid. */
  718.  
  719.               h = P[Y] + k * D[Y];
  720.  
  721.               if ((h >= Prism->Height1) && (h <= Prism->Height2))
  722.               {
  723.                 Intersection[i].t   = SPLINE_HIT;
  724.                 Intersection[i].n   = j;
  725.                 Intersection[i].w   = w;
  726.                 Intersection[i++].d = k / len;
  727.               }
  728.             }
  729.           }
  730.         }
  731.       }
  732.  
  733.       break;
  734.  
  735.       default:
  736.  
  737.         Error("Unknown sweep type in intersect_prism().\n");
  738.   }
  739.  
  740.   if (i)
  741.   {
  742.     Increase_Counter(stats[Ray_Prism_Tests_Succeeded]);
  743.   }
  744.  
  745.   return(i);
  746. }
  747.  
  748.  
  749.  
  750. /*****************************************************************************
  751. *
  752. * FUNCTION
  753. *
  754. *   Inside_Prism
  755. *
  756. * INPUT
  757. *
  758. *   IPoint - Intersection point
  759. *   Object - Object
  760. *   
  761. * OUTPUT
  762. *   
  763. * RETURNS
  764. *
  765. *   int - TRUE if inside
  766. *   
  767. * AUTHOR
  768. *
  769. *   Dieter Bayer
  770. *   
  771. * DESCRIPTION
  772. *
  773. *   Test if point lies inside a prism.
  774. *
  775. * CHANGES
  776. *
  777. *   May 1994 : Creation.
  778. *
  779. ******************************************************************************/
  780.  
  781. static int Inside_Prism(IPoint, Object)
  782. VECTOR IPoint;
  783. OBJECT *Object;
  784. {
  785.   VECTOR P;
  786.   PRISM *Prism = (PRISM *)Object;
  787.   DBL rdiv;
  788.   
  789.   /* Transform the point into the prism space. */
  790.  
  791.   MInvTransPoint(P, IPoint, Prism->Trans);
  792.  
  793.   if ((P[Y] >= Prism->Height1) && (P[Y] < Prism->Height2))
  794.   {
  795.     if (Prism->Sweep_Type == CONIC_SWEEP)
  796.     {
  797.       /* Scale x and z coordinate. */
  798.  
  799.       if (P[Y] > 0.0)
  800.       {
  801.         rdiv = (1.0 / P[Y]);
  802.         P[X] *= rdiv;
  803.         P[Z] *= rdiv;
  804. /*
  805.         P[X] /= P[Y];
  806.         P[Z] /= P[Y];
  807. */
  808.       }
  809.       else
  810.       {
  811.         P[X] = P[Z] = HUGE_VAL;
  812.       }
  813.     }
  814.  
  815.     if (in_curve(Prism, P[X], P[Z]))
  816.     {
  817.       return(!Test_Flag(Prism, INVERTED_FLAG));
  818.     }
  819.   }
  820.  
  821.   return(Test_Flag(Prism, INVERTED_FLAG));
  822. }
  823.  
  824.  
  825.  
  826. /*****************************************************************************
  827. *
  828. * FUNCTION
  829. *
  830. *   Prism_Normal
  831. *
  832. * INPUT
  833. *
  834. *   Result - Normal vector
  835. *   Object - Object
  836. *   Inter  - Intersection found
  837. *   
  838. * OUTPUT
  839. *
  840. *   Result
  841. *   
  842. * RETURNS
  843. *   
  844. * AUTHOR
  845. *
  846. *   Dieter Bayer
  847. *   
  848. * DESCRIPTION
  849. *
  850. *   Calculate the normal of the prism in a given point.
  851. *
  852. * CHANGES
  853. *
  854. *   May 1994 : Creation.
  855. *
  856. ******************************************************************************/
  857.  
  858. static void Prism_Normal(Result, Object, Inter)
  859. OBJECT *Object;
  860. VECTOR Result;
  861. INTERSECTION *Inter;
  862. {
  863.   DBL r;
  864.   VECTOR P;
  865.   PRISM_SPLINE_ENTRY Entry;
  866.   PRISM *Prism = (PRISM *)Object;
  867.   VECTOR N;
  868.  
  869.   Make_Vector(N, 0.0, 1.0, 0.0);
  870.  
  871.   if (Inter->i1 == SPLINE_HIT)
  872.   {
  873.     Entry = Prism->Spline->Entry[Inter->i2];
  874.  
  875.     switch (Prism->Sweep_Type)
  876.     {
  877.       case LINEAR_SWEEP:
  878.  
  879.         N[X] =   Inter->d1 * (3.0 * Entry.A[Y] * Inter->d1 + 2.0 * Entry.B[Y]) + Entry.C[Y];
  880.         N[Y] =   0.0;
  881.         N[Z] = -(Inter->d1 * (3.0 * Entry.A[X] * Inter->d1 + 2.0 * Entry.B[X]) + Entry.C[X]);
  882.  
  883.         break;
  884.  
  885.       case CONIC_SWEEP:
  886.  
  887.         /* Transform the point into the prism space. */
  888.  
  889.         MInvTransPoint(P, Inter->IPoint, Prism->Trans);
  890.  
  891.         if (P[Y] > 0.0)
  892.         {
  893.           r = sqrt(P[X] * P[X] + P[Z] * P[Z]);
  894.  
  895.           N[X] =  Inter->d1 * (3.0 * Entry.A[Y] * Inter->d1 + 2.0 * Entry.B[Y]) + Entry.C[Y];
  896.           N[Y] = -r / P[Y];
  897.           N[Z] = -(Inter->d1 * (3.0 * Entry.A[X] * Inter->d1 + 2.0 * Entry.B[X]) + Entry.C[X]);
  898.         }
  899.  
  900.         break;
  901.  
  902.       default:
  903.  
  904.         Error("Unknown sweep type in Prism_Normal().\n");
  905.     }
  906.   }
  907.  
  908.   /* Transform the normalt out of the prism space. */
  909.  
  910.   MTransNormal(Result, N, Prism->Trans);
  911.  
  912.   VNormalize(Result, Result);
  913. }
  914.  
  915.  
  916.  
  917. /*****************************************************************************
  918. *
  919. * FUNCTION
  920. *
  921. *   Translate_Prism
  922. *
  923. * INPUT
  924. *
  925. *   Object - Object
  926. *   Vector - Translation vector
  927. *   
  928. * OUTPUT
  929. *
  930. *   Object
  931. *   
  932. * RETURNS
  933. *   
  934. * AUTHOR
  935. *
  936. *   Dieter Bayer
  937. *   
  938. * DESCRIPTION
  939. *
  940. *   Translate a prism.
  941. *
  942. * CHANGES
  943. *
  944. *   May 1994 : Creation.
  945. *
  946. ******************************************************************************/
  947.  
  948. static void Translate_Prism(Object, Vector, Trans)
  949. OBJECT *Object;
  950. VECTOR Vector;
  951. TRANSFORM *Trans;
  952. {
  953.   Transform_Prism(Object, Trans);
  954. }
  955.  
  956.  
  957.  
  958. /*****************************************************************************
  959. *
  960. * FUNCTION
  961. *
  962. *   Rotate_Prism
  963. *
  964. * INPUT
  965. *
  966. *   Object - Object
  967. *   Vector - Rotation vector
  968. *   
  969. * OUTPUT
  970. *
  971. *   Object
  972. *   
  973. * RETURNS
  974. *   
  975. * AUTHOR
  976. *
  977. *   Dieter Bayer
  978. *   
  979. * DESCRIPTION
  980. *
  981. *   Rotate a prism.
  982. *
  983. * CHANGES
  984. *
  985. *   May 1994 : Creation.
  986. *
  987. ******************************************************************************/
  988.  
  989. static void Rotate_Prism(Object, Vector, Trans)
  990. OBJECT *Object;
  991. VECTOR Vector;
  992. TRANSFORM *Trans;
  993. {
  994.   Transform_Prism(Object, Trans);
  995. }
  996.  
  997.  
  998.  
  999. /*****************************************************************************
  1000. *
  1001. * FUNCTION
  1002. *
  1003. *   Scale_Prism
  1004. *
  1005. * INPUT
  1006. *
  1007. *   Object - Object
  1008. *   Vector - Scaling vector
  1009. *   
  1010. * OUTPUT
  1011. *
  1012. *   Object
  1013. *   
  1014. * RETURNS
  1015. *   
  1016. * AUTHOR
  1017. *
  1018. *   Dieter Bayer
  1019. *
  1020. * DESCRIPTION
  1021. *
  1022. *   Scale a prism.
  1023. *
  1024. * CHANGES
  1025. *
  1026. *   May 1994 : Creation.
  1027. *
  1028. ******************************************************************************/
  1029.  
  1030. static void Scale_Prism(Object, Vector, Trans)
  1031. OBJECT *Object;
  1032. VECTOR Vector;
  1033. TRANSFORM *Trans;
  1034. {
  1035.   Transform_Prism(Object, Trans);
  1036. }
  1037.  
  1038.  
  1039.  
  1040. /*****************************************************************************
  1041. *
  1042. * FUNCTION
  1043. *
  1044. *   Transform_Prism
  1045. *
  1046. * INPUT
  1047. *
  1048. *   Object - Object
  1049. *   Trans  - Transformation to apply
  1050. *   
  1051. * OUTPUT
  1052. *
  1053. *   Object
  1054. *   
  1055. * RETURNS
  1056. *   
  1057. * AUTHOR
  1058. *
  1059. *   Dieter Bayer
  1060. *   
  1061. * DESCRIPTION
  1062. *
  1063. *   Transform a prism and recalculate its bounding box.
  1064. *
  1065. * CHANGES
  1066. *
  1067. *   May 1994 : Creation.
  1068. *
  1069. ******************************************************************************/
  1070.  
  1071. static void Transform_Prism(Object, Trans)
  1072. OBJECT *Object;
  1073. TRANSFORM *Trans;
  1074. {
  1075.   Compose_Transforms(((PRISM *)Object)->Trans, Trans);
  1076.  
  1077.   Compute_Prism_BBox((PRISM *)Object);
  1078. }
  1079.  
  1080.  
  1081.  
  1082. /*****************************************************************************
  1083. *
  1084. * FUNCTION
  1085. *
  1086. *   Invert_Prism
  1087. *
  1088. * INPUT
  1089. *
  1090. *   Object - Object
  1091. *   
  1092. * OUTPUT
  1093. *
  1094. *   Object
  1095. *   
  1096. * RETURNS
  1097. *   
  1098. * AUTHOR
  1099. *
  1100. *   Dieter Bayer
  1101. *   
  1102. * DESCRIPTION
  1103. *
  1104. *   Invert a prism.
  1105. *
  1106. * CHANGES
  1107. *
  1108. *   May 1994 : Creation.
  1109. *
  1110. ******************************************************************************/
  1111.  
  1112. static void Invert_Prism(Object)
  1113. OBJECT *Object;
  1114. {
  1115.   Invert_Flag(Object, INVERTED_FLAG);
  1116. }
  1117.  
  1118.  
  1119.  
  1120. /*****************************************************************************
  1121. *
  1122. * FUNCTION
  1123. *
  1124. *   Create_Prism
  1125. *
  1126. * INPUT
  1127. *   
  1128. * OUTPUT
  1129. *   
  1130. * RETURNS
  1131. *
  1132. *   PRISM * - new prism
  1133. *   
  1134. * AUTHOR
  1135. *
  1136. *   Dieter Bayer
  1137. *   
  1138. * DESCRIPTION
  1139. *
  1140. *   Create a new prism.
  1141. *
  1142. * CHANGES
  1143. *
  1144. *   May 1994 : Creation.
  1145. *
  1146. ******************************************************************************/
  1147.  
  1148. PRISM *Create_Prism()
  1149. {
  1150.   PRISM *New;
  1151.  
  1152.   New = (PRISM *)POV_MALLOC(sizeof(PRISM), "prism");
  1153.  
  1154.   INIT_OBJECT_FIELDS(New,PRISM_OBJECT,&Prism_Methods)
  1155.  
  1156.   New->Trans = Create_Transform();
  1157.  
  1158.   New->x1      =
  1159.   New->x2      =
  1160.   New->y1      =
  1161.   New->y2      =
  1162.   New->Height1 =
  1163.   New->Height2 = 0.0;
  1164.  
  1165.   New->Number = 0;
  1166.  
  1167.   New->Spline_Type = LINEAR_SPLINE;
  1168.   New->Sweep_Type  = LINEAR_SWEEP;
  1169.  
  1170.   New->Spline = NULL;
  1171.  
  1172.   Set_Flag(New, CLOSED_FLAG);
  1173.  
  1174.   return(New);
  1175. }
  1176.  
  1177.  
  1178.  
  1179. /*****************************************************************************
  1180. *
  1181. * FUNCTION
  1182. *
  1183. *   Copy_Prism
  1184. *
  1185. * INPUT
  1186. *
  1187. *   Object - Object
  1188. *   
  1189. * OUTPUT
  1190. *   
  1191. * RETURNS
  1192. *
  1193. *   void * - New prism
  1194. *   
  1195. * AUTHOR
  1196. *
  1197. *   Dieter Bayer
  1198. *   
  1199. * DESCRIPTION
  1200. *
  1201. *   Copy a prism structure.
  1202. *
  1203. *   NOTE: The splines are not copied, only the number of references is
  1204. *         counted, so that Destray_Prism() knows if they can be destroyed.
  1205. *
  1206. * CHANGES
  1207. *
  1208. *   May 1994 : Creation.
  1209. *
  1210. *   Sep 1994 : fixed memory leakage [DB]
  1211. *
  1212. ******************************************************************************/
  1213.  
  1214. static void *Copy_Prism(Object)
  1215. OBJECT *Object;
  1216. {
  1217.   PRISM *New, *Prism = (PRISM *)Object;
  1218.  
  1219.   New = Create_Prism();
  1220.  
  1221.   /* Get rid of the transformation created in Create_Prism(). */
  1222.  
  1223.   Destroy_Transform(New->Trans);
  1224.  
  1225.   /* Copy prism. */
  1226.  
  1227.   *New = *Prism;
  1228.  
  1229.   New->Trans = Copy_Transform(((PRISM *)Object)->Trans);
  1230.  
  1231.   New->Spline->References++;
  1232.  
  1233.   return(New);
  1234. }
  1235.  
  1236.  
  1237.  
  1238. /*****************************************************************************
  1239. *
  1240. * FUNCTION
  1241. *
  1242. *   Destroy_Prism
  1243. *
  1244. * INPUT
  1245. *
  1246. *   Object - Object
  1247. *   
  1248. * OUTPUT
  1249. *
  1250. *   Object
  1251. *   
  1252. * RETURNS
  1253. *   
  1254. * AUTHOR
  1255. *
  1256. *   Dieter Bayer
  1257. *   
  1258. * DESCRIPTION
  1259. *
  1260. *   Destroy a prism.
  1261. *
  1262. *   NOTE: The splines are destroyed if they are no longer used by any copy.
  1263. *
  1264. * CHANGES
  1265. *
  1266. *   May 1994 : Creation.
  1267. *
  1268. ******************************************************************************/
  1269.  
  1270. static void Destroy_Prism (Object)
  1271. OBJECT *Object;
  1272. {
  1273.   PRISM *Prism = (PRISM *)Object;
  1274.  
  1275.   Destroy_Transform(Prism->Trans);
  1276.  
  1277.   if (--(Prism->Spline->References) == 0)
  1278.   {
  1279.     POV_FREE (Prism->Spline->Entry);
  1280.  
  1281.     POV_FREE (Prism->Spline);
  1282.   }
  1283.  
  1284.   POV_FREE (Object);
  1285. }
  1286.  
  1287.  
  1288.  
  1289. /*****************************************************************************
  1290. *
  1291. * FUNCTION
  1292. *
  1293. *   Compute_Prism_BBox
  1294. *
  1295. * INPUT
  1296. *
  1297. *   Prism - Prism
  1298. *   
  1299. * OUTPUT
  1300. *
  1301. *   Prism
  1302. *   
  1303. * RETURNS
  1304. *   
  1305. * AUTHOR
  1306. *
  1307. *   Dieter Bayer
  1308. *   
  1309. * DESCRIPTION
  1310. *
  1311. *   Calculate the bounding box of a prism.
  1312. *
  1313. * CHANGES
  1314. *
  1315. *   May 1994 : Creation.
  1316. *
  1317. ******************************************************************************/
  1318.  
  1319. void Compute_Prism_BBox(Prism)
  1320. PRISM *Prism;
  1321. {
  1322.   Make_BBox(Prism->BBox, Prism->x1, Prism->Height1, Prism->y1,
  1323.     Prism->x2 - Prism->x1, Prism->Height2 - Prism->Height1, Prism->y2 - Prism->y1);
  1324.  
  1325.   Recompute_BBox(&Prism->BBox, Prism->Trans);
  1326. }
  1327.  
  1328.  
  1329.  
  1330. /*****************************************************************************
  1331. *
  1332. * FUNCTION
  1333. *
  1334. *   in_curve
  1335. *
  1336. * INPUT
  1337. *
  1338. *   Prism - Prism to test
  1339. *   u, v  - Coordinates
  1340. *   
  1341. * OUTPUT
  1342. *   
  1343. * RETURNS
  1344. *
  1345. *   int - TRUE if inside
  1346. *   
  1347. * AUTHOR
  1348. *
  1349. *   Dieter Bayer, June 1994
  1350. *   
  1351. * DESCRIPTION
  1352. *
  1353. *   Test if a 2d point lies inside a prism's spline curve.
  1354. *
  1355. *   We travel from the current point in positive u direction
  1356. *   and count the number of crossings with the curve.
  1357. *
  1358. * CHANGES
  1359. *
  1360. *   May 1994 : Creation.
  1361. *
  1362. ******************************************************************************/
  1363.  
  1364. static int in_curve(Prism, u, v)
  1365. PRISM *Prism;
  1366. DBL u, v;
  1367. {
  1368.   int i, n, NC;
  1369.   DBL k, w;
  1370.   DBL x[4], y[3];
  1371.   PRISM_SPLINE_ENTRY Entry;
  1372.  
  1373.   NC = 0;
  1374.  
  1375.   /* First test overall bounding rectangle. */
  1376.   
  1377.   if ((u >= Prism->x1) && (u <= Prism->x2) &&
  1378.       (v >= Prism->y1) && (v <= Prism->y2))
  1379.   {
  1380.     for (i = 0; i < Prism->Number; i++)
  1381.     {
  1382.       Entry = Prism->Spline->Entry[i];
  1383.  
  1384.       /* Test if current segment can be hit. */
  1385.  
  1386.       if ((v >= Entry.y1) && (v <= Entry.y2) && (u <= Entry.x2))
  1387.       {
  1388.         x[0] = Entry.A[Y];
  1389.         x[1] = Entry.B[Y];
  1390.         x[2] = Entry.C[Y];
  1391.         x[3] = Entry.D[Y] - v;
  1392.  
  1393.         n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
  1394.  
  1395.         while (n--)
  1396.         {
  1397.           w = y[n];
  1398.  
  1399.           if ((w >= 0.0) && (w <= 1.0))
  1400.           {
  1401.             k  = w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X] - u;
  1402.  
  1403.             if (k >= 0.0)
  1404.             {
  1405.               NC++;
  1406.             }
  1407.           }
  1408.         }
  1409.       }
  1410.     }
  1411.   }
  1412.  
  1413.   return(NC & 1);
  1414. }
  1415.  
  1416.  
  1417.  
  1418. /*****************************************************************************
  1419. *
  1420. * FUNCTION
  1421. *
  1422. *   test_rectangle
  1423. *
  1424. * INPUT
  1425. *   
  1426. * OUTPUT
  1427. *   
  1428. * RETURNS
  1429. *   
  1430. * AUTHOR
  1431. *
  1432. *   Dieter Bayer, July 1994
  1433. *   
  1434. * DESCRIPTION
  1435. *
  1436. *   Test if the 2d ray (= P + t * D) intersects a rectangle.
  1437. *
  1438. * CHANGES
  1439. *
  1440. *   May 1994 : Creation.
  1441. *
  1442. ******************************************************************************/
  1443.  
  1444. static int test_rectangle(P, D, x1, z1, x2, z2)
  1445. VECTOR P, D;
  1446. DBL x1, z1, x2, z2;
  1447. {
  1448.   DBL dmin, dmax, tmin, tmax;
  1449.   DBL rdivx;
  1450.   DBL rdivz;
  1451.   
  1452.   if (fabs(D[X]) > EPSILON)
  1453.   {
  1454.   rdivx = (1.0 / D[X]);
  1455.  
  1456.     if (D[X] > 0.0)
  1457.     {
  1458.       dmin = (x1 - P[X]) * rdivx;
  1459.       dmax = (x2 - P[X]) * rdivx;
  1460. /*
  1461.       dmin = (x1 - P[X]) / D[X];
  1462.       dmax = (x2 - P[X]) / D[X];
  1463. */
  1464.       if (dmax < EPSILON)
  1465.       {
  1466.         return(FALSE);
  1467.       }
  1468.     }
  1469.     else
  1470.     {
  1471.       dmax = (x1 - P[X]) * rdivx;
  1472. /*
  1473.       dmax = (x1 - P[X]) / D[X];
  1474. */
  1475.       if (dmax < EPSILON)
  1476.       {
  1477.         return(FALSE);
  1478.       }
  1479.       dmin = (x2 - P[X]) * rdivx;
  1480. /*
  1481.       dmin = (x2 - P[X]) / D[X];
  1482. */
  1483.     }
  1484.  
  1485.     if (dmin > dmax)
  1486.     {
  1487.       return(FALSE);
  1488.     }
  1489.   }
  1490.   else
  1491.   {
  1492.     if ((P[X] < x1) || (P[X] > x2))
  1493.     {
  1494.       return(FALSE);
  1495.     }
  1496.     else
  1497.     {
  1498.       dmin = -BOUND_HUGE;
  1499.       dmax =  BOUND_HUGE;
  1500.     }
  1501.   }
  1502.  
  1503.   if (fabs(D[Z]) > EPSILON)
  1504.   {
  1505.     rdivz = (1.0 / D[Z]);
  1506.     if (D[Z] > 0.0)
  1507.     {
  1508.       tmin = (z1 - P[Z]) * rdivz;
  1509.       tmax = (z2 - P[Z]) * rdivz;
  1510. /*
  1511.       tmin = (z1 - P[Z]) / D[Z];
  1512.       tmax = (z2 - P[Z]) / D[Z];
  1513. */
  1514.     }
  1515.     else
  1516.     {
  1517.       tmax = (z1 - P[Z]) * rdivz;
  1518.       tmin = (z2 - P[Z]) * rdivz;
  1519. /*
  1520.       tmax = (z1 - P[Z]) / D[Z];
  1521.       tmin = (z2 - P[Z]) / D[Z];
  1522. */
  1523.     }
  1524.  
  1525.     if (tmax < dmax)
  1526.     {
  1527.       if (tmax < EPSILON)
  1528.       {
  1529.         return(FALSE);
  1530.       }
  1531.  
  1532.       if (tmin > dmin)
  1533.       {
  1534.         if (tmin > tmax)
  1535.         {
  1536.           return(FALSE);
  1537.         }
  1538.  
  1539.         dmin = tmin;
  1540.       }
  1541.       else
  1542.       {
  1543.         if (dmin > tmax)
  1544.         {
  1545.           return(FALSE);
  1546.         }
  1547.       }
  1548.  
  1549.       /*dmax = tmax; */ /*not needed CEY[1/95]*/
  1550.     }
  1551.     else
  1552.     {
  1553.       if (tmin > dmin)
  1554.       {
  1555.         if (tmin > dmax)
  1556.         {
  1557.           return(FALSE);
  1558.         }
  1559.  
  1560.         /* dmin = tmin; */  /*not needed CEY[1/95]*/
  1561.       }
  1562.     }
  1563.   }
  1564.   else
  1565.   {
  1566.     if ((P[Z] < z1) || (P[Z] > z2))
  1567.     {
  1568.       return(FALSE);
  1569.     }
  1570.   }
  1571.  
  1572.   return(TRUE);
  1573. }
  1574.  
  1575.  
  1576.  
  1577. /*****************************************************************************
  1578. *
  1579. * FUNCTION
  1580. *
  1581. *   Compute_Prism
  1582. *
  1583. * INPUT
  1584. *
  1585. *   Prism - Prism
  1586. *   P     - Points defining prism
  1587. *   
  1588. * OUTPUT
  1589. *
  1590. *   Prism
  1591. *
  1592. * RETURNS
  1593. *
  1594. * AUTHOR
  1595. *
  1596. *   Dieter Bayer, June 1994
  1597. *
  1598. * DESCRIPTION
  1599. *
  1600. *   Calculate the spline segments of a prism from a set of points.
  1601. *
  1602. * CHANGES
  1603. *
  1604. *   May 1994 : Creation.
  1605. *
  1606. ******************************************************************************/
  1607.  
  1608. void Compute_Prism(Prism, P)
  1609. PRISM *Prism;
  1610. UV_VECT *P;
  1611. {
  1612.   int i, n, number_of_splines;
  1613.   int i1, i2, i3;
  1614.  
  1615.   DBL x[4], xmin, xmax;
  1616.   DBL y[4], ymin, ymax;
  1617.   DBL c[3], r[2];
  1618.  
  1619.   UV_VECT A, B, C, D, First;
  1620.  
  1621.   /* Allocate Object->Number segments. */
  1622.  
  1623.   if (Prism->Spline == NULL)
  1624.   {
  1625.     Prism->Spline = (PRISM_SPLINE *)POV_MALLOC(sizeof(PRISM_SPLINE), "spline segments of prism");
  1626.  
  1627.     Prism->Spline->References = 1;
  1628.  
  1629.     Prism->Spline->Entry = (PRISM_SPLINE_ENTRY *)POV_MALLOC(Prism->Number*sizeof(PRISM_SPLINE_ENTRY), "spline segments of prism");
  1630.   }
  1631.   else
  1632.   {
  1633.     /* This should never happen! */
  1634.  
  1635.     Error("Prism segments are already defined.\n");
  1636.   }
  1637.  
  1638.   /***************************************************************************
  1639.   * Calculate segments.
  1640.   ****************************************************************************/
  1641.  
  1642.   /* We want to know the size of the overall bounding rectangle. */
  1643.  
  1644.   xmax = ymax = -BOUND_HUGE;
  1645.   xmin = ymin =  BOUND_HUGE;
  1646.  
  1647.   /* Get first segment point. */
  1648.  
  1649.   switch (Prism->Spline_Type)
  1650.   {
  1651.     case LINEAR_SPLINE:
  1652.       Assign_UV_Vect(First, P[0]);
  1653.  
  1654.       break;
  1655.  
  1656.     case QUADRATIC_SPLINE:
  1657.     case CUBIC_SPLINE:
  1658.  
  1659.       Assign_UV_Vect(First, P[1]);
  1660.  
  1661.       break;
  1662.   }
  1663.  
  1664.   for (i = number_of_splines = 0; i < Prism->Number-1; i++)
  1665.   {
  1666.     /* Get indices of previous and next two points. */
  1667.  
  1668.     i1 = i + 1;
  1669.     i2 = i + 2;
  1670.     i3 = i + 3;
  1671.  
  1672.     switch (Prism->Spline_Type)
  1673.     {
  1674.       /*************************************************************************
  1675.       * Linear spline (nothing more than a simple polygon).
  1676.       **************************************************************************/
  1677.  
  1678.       case LINEAR_SPLINE :
  1679.  
  1680.         if (i1 >= Prism->Number)
  1681.         {
  1682.           Error("Too few points in prism. Prism not closed? Control points missing?\n");
  1683.         }
  1684.  
  1685.         /* Use linear interpolation. */
  1686.  
  1687.         A[X] =  0.0;
  1688.         B[X] =  0.0;
  1689.         C[X] = -1.0 * P[i][X] + 1.0 * P[i1][X];
  1690.         D[X] =  1.0 * P[i][X];
  1691.  
  1692.         A[Y] =  0.0;
  1693.         B[Y] =  0.0;
  1694.         C[Y] = -1.0 * P[i][Y] + 1.0 * P[i1][Y];
  1695.         D[Y] =  1.0 * P[i][Y];
  1696.  
  1697.         x[0] = x[2] = P[i][X];
  1698.         x[1] = x[3] = P[i1][X];
  1699.  
  1700.         y[0] = y[2] = P[i][Y];
  1701.         y[1] = y[3] = P[i1][Y];
  1702.  
  1703.         break;
  1704.  
  1705.       /*************************************************************************
  1706.       * Quadratic spline.
  1707.       **************************************************************************/
  1708.  
  1709.       case QUADRATIC_SPLINE :
  1710.  
  1711.         if (i2 >= Prism->Number)
  1712.         {
  1713.           Error("Too few points in prism.\n");
  1714.         }
  1715.  
  1716.         /* Use quadratic interpolation. */
  1717.  
  1718.         A[X] =  0.0;
  1719.         B[X] =  0.5 * P[i][X] - 1.0 * P[i1][X] + 0.5 * P[i2][X];
  1720.         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
  1721.         D[X] =                  1.0 * P[i1][X];
  1722.  
  1723.         A[Y] =  0.0;
  1724.         B[Y] =  0.5 * P[i][Y] - 1.0 * P[i1][Y] + 0.5 * P[i2][Y];
  1725.         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
  1726.         D[Y] =                  1.0 * P[i1][Y];
  1727.  
  1728.         x[0] = x[2] = P[i1][X];
  1729.         x[1] = x[3] = P[i2][X];
  1730.  
  1731.         y[0] = y[2] = P[i1][Y];
  1732.         y[1] = y[3] = P[i2][Y];
  1733.  
  1734.         break;
  1735.  
  1736.       /*************************************************************************
  1737.       * Cubic spline.
  1738.       **************************************************************************/
  1739.  
  1740.       case CUBIC_SPLINE :
  1741.  
  1742.         if (i3 >= Prism->Number)
  1743.         {
  1744.           Error("Too few points in prism.\n");
  1745.         }
  1746.  
  1747.         /* Use cubic interpolation. */
  1748.  
  1749.         A[X] = -0.5 * P[i][X] + 1.5 * P[i1][X] - 1.5 * P[i2][X] + 0.5 * P[i3][X];
  1750.         B[X] =        P[i][X] - 2.5 * P[i1][X] + 2.0 * P[i2][X] - 0.5 * P[i3][X];
  1751.         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
  1752.         D[X] =                        P[i1][X];
  1753.  
  1754.         A[Y] = -0.5 * P[i][Y] + 1.5 * P[i1][Y] - 1.5 * P[i2][Y] + 0.5 * P[i3][Y];
  1755.         B[Y] =        P[i][Y] - 2.5 * P[i1][Y] + 2.0 * P[i2][Y] - 0.5 * P[i3][Y];
  1756.         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
  1757.         D[Y] =                        P[i1][Y];
  1758.  
  1759.         x[0] = x[2] = P[i1][X];
  1760.         x[1] = x[3] = P[i2][X];
  1761.  
  1762.         y[0] = y[2] = P[i1][Y];
  1763.         y[1] = y[3] = P[i2][Y];
  1764.  
  1765.         break;
  1766.  
  1767.       default:
  1768.  
  1769.         Error("Unknown spline type in Compute_Prism().\n");
  1770.     }
  1771.  
  1772.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].A, A);
  1773.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].B, B);
  1774.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].C, C);
  1775.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].D, D);
  1776.  
  1777.     /* Get maximum coordinates in current segment. */
  1778.  
  1779.     c[0] = 3.0 * A[X];
  1780.     c[1] = 2.0 * B[X];
  1781.     c[2] = C[X];
  1782.  
  1783.     n = Solve_Polynomial(2, c, r, FALSE, 0.0);
  1784.  
  1785.     while (n--)
  1786.     {
  1787.       if ((r[n] >= 0.0) && (r[n] <= 1.0))
  1788.       {
  1789.         x[n] = r[n] * (r[n] * (r[n] * A[X] + B[X]) + C[X]) + D[X];
  1790.       }
  1791.     }
  1792.  
  1793.     c[0] = 3.0 * A[Y];
  1794.     c[1] = 2.0 * B[Y];
  1795.     c[2] = C[Y];
  1796.  
  1797.     n = Solve_Polynomial(2, c, r, FALSE, 0.0);
  1798.  
  1799.     while (n--)
  1800.     {
  1801.       if ((r[n] >= 0.0) && (r[n] <= 1.0))
  1802.       {
  1803.         y[n] = r[n] * (r[n] * (r[n] * A[Y] + B[Y]) + C[Y]) + D[Y];
  1804.       }
  1805.     }
  1806.  
  1807.     /* Set current segment's bounding rectangle. */
  1808.  
  1809.     Prism->Spline->Entry[number_of_splines].x1 = min(min(x[0], x[1]), min(x[2], x[3]));
  1810.     Prism->Spline->Entry[number_of_splines].x2 = max(max(x[0], x[1]), max(x[2], x[3]));
  1811.  
  1812.     Prism->Spline->Entry[number_of_splines].y1 = min(min(y[0], y[1]), min(y[2], y[3]));
  1813.     Prism->Spline->Entry[number_of_splines].y2 = max(max(y[0], y[1]), max(y[2], y[3]));
  1814.  
  1815.     /* Keep track of overall bounding rectangle. */
  1816.  
  1817.     xmin = min(xmin, Prism->Spline->Entry[number_of_splines].x1);
  1818.     xmax = max(xmax, Prism->Spline->Entry[number_of_splines].x2);
  1819.  
  1820.     ymin = min(ymin, Prism->Spline->Entry[number_of_splines].y1);
  1821.     ymax = max(ymax, Prism->Spline->Entry[number_of_splines].y2);
  1822.  
  1823.     number_of_splines++;
  1824.  
  1825.     /* Check if end of sub-prism is reached. */
  1826.  
  1827.     switch (Prism->Spline_Type)
  1828.     {
  1829.       case LINEAR_SPLINE:
  1830.  
  1831.         if ((fabs(P[i1][X] - First[X]) < EPSILON) &&
  1832.             (fabs(P[i1][Y] - First[Y]) < EPSILON))
  1833.         {
  1834.           i++;
  1835.  
  1836.           if (i+1 < Prism->Number)
  1837.           {
  1838.             Assign_UV_Vect(First, P[i+1]);
  1839.           }
  1840.         }
  1841.  
  1842.         break;
  1843.  
  1844.       case QUADRATIC_SPLINE:
  1845.  
  1846.         if ((fabs(P[i2][X] - First[X]) < EPSILON) &&
  1847.             (fabs(P[i2][Y] - First[Y]) < EPSILON))
  1848.         {
  1849.           i += 2;
  1850.  
  1851.           if (i+2 < Prism->Number)
  1852.           {
  1853.             Assign_UV_Vect(First, P[i+2]);
  1854.           }
  1855.         }
  1856.  
  1857.         break;
  1858.  
  1859.       case CUBIC_SPLINE:
  1860.  
  1861.         if ((fabs(P[i2][X] - First[X]) < EPSILON) &&
  1862.             (fabs(P[i2][Y] - First[Y]) < EPSILON))
  1863.         {
  1864.           i += 3;
  1865.  
  1866.           if (i+2 < Prism->Number)
  1867.           {
  1868.             Assign_UV_Vect(First, P[i+2]);
  1869.           }
  1870.         }
  1871.  
  1872.         break;
  1873.  
  1874.     }
  1875.   }
  1876.  
  1877.   Prism->Number = number_of_splines;
  1878.  
  1879. /*
  1880.   Prism->Spline->Entry = (PRISM_SPLINE_ENTRY *)POV_REALLOC(Prism->Spline->Entry,
  1881.     Prism->Number*sizeof(PRISM_SPLINE_ENTRY), "spline segments of prism");
  1882. */
  1883.  
  1884.   /* Set overall bounding rectangle. */
  1885.  
  1886.   Prism->x1 = xmin;
  1887.   Prism->x2 = xmax;
  1888.  
  1889.   Prism->y1 = ymin;
  1890.   Prism->y2 = ymax;
  1891.  
  1892.   if (Prism->Sweep_Type == CONIC_SWEEP)
  1893.   {
  1894.     /* Recalculate bounding rectangles. */
  1895.  
  1896.     for (i = 0; i < Prism->Number; i++)
  1897.     {
  1898.       x[0] = Prism->Spline->Entry[i].x1 * Prism->Height1;
  1899.       x[1] = Prism->Spline->Entry[i].x1 * Prism->Height2;
  1900.       x[2] = Prism->Spline->Entry[i].x2 * Prism->Height1;
  1901.       x[3] = Prism->Spline->Entry[i].x2 * Prism->Height2;
  1902.  
  1903.       xmin = min(min(x[0], x[1]), min(x[2], x[3]));
  1904.       xmax = max(max(x[0], x[1]), max(x[2], x[3]));
  1905.  
  1906.       Prism->Spline->Entry[i].x1 = xmin;
  1907.       Prism->Spline->Entry[i].x2 = xmax;
  1908.  
  1909.       y[0] = Prism->Spline->Entry[i].y1 * Prism->Height1;
  1910.       y[1] = Prism->Spline->Entry[i].y1 * Prism->Height2;
  1911.       y[2] = Prism->Spline->Entry[i].y2 * Prism->Height1;
  1912.       y[3] = Prism->Spline->Entry[i].y2 * Prism->Height2;
  1913.  
  1914.       ymin = min(min(y[0], y[1]), min(y[2], y[3]));
  1915.       ymax = max(max(y[0], y[1]), max(y[2], y[3]));
  1916.  
  1917.       Prism->Spline->Entry[i].y1 = ymin;
  1918.       Prism->Spline->Entry[i].y2 = ymax;
  1919.     }
  1920.  
  1921.     /* Recalculate overall bounding rectangle. */
  1922.  
  1923.     x[0] = Prism->x1 * Prism->Height1;
  1924.     x[1] = Prism->x1 * Prism->Height2;
  1925.     x[2] = Prism->x2 * Prism->Height1;
  1926.     x[3] = Prism->x2 * Prism->Height2;
  1927.  
  1928.     xmin = min(min(x[0], x[1]), min(x[2], x[3]));
  1929.     xmax = max(max(x[0], x[1]), max(x[2], x[3]));
  1930.  
  1931.     Prism->x1 = xmin;
  1932.     Prism->x2 = xmax;
  1933.  
  1934.     y[0] = Prism->y1 * Prism->Height1;
  1935.     y[1] = Prism->y1 * Prism->Height2;
  1936.     y[2] = Prism->y2 * Prism->Height1;
  1937.     y[3] = Prism->y2 * Prism->Height2;
  1938.  
  1939.     ymin = min(min(y[0], y[1]), min(y[2], y[3]));
  1940.     ymax = max(max(y[0], y[1]), max(y[2], y[3]));
  1941.  
  1942.     Prism->y1 = ymin;
  1943.     Prism->y2 = ymax;
  1944.   }
  1945. }
  1946.  
  1947.  
  1948.  
  1949.